Distance  Functions  and  Geodesics  on  Points  Clouds 


Facundo  Memoli*  Guillermo  Sapiro1 


Abstract 

An  algorithm  for  computing  intrinsic  distance  functions  and  geodesics  on  sub-manifolds  of  Md  given 
by  point  clouds  is  introduced  in  this  paper.  The  basic  idea  is  that,  as  shown  in  this  paper,  intrinsic  distance 
functions  and  geodesics  on  general  co-dimension  sub-manifolds  of  Md  can  be  accurately  approximated 
by  the  extrinsic  Euclidean  ones  computed  in  a  thin  offset  band  surrounding  the  manifold.  This  permits 
the  use  of  computationally  optimal  algorithms  for  computing  distance  functions  in  Cartesian  grids.  We 
then  use  these  algorithms,  modified  to  deal  with  spaces  with  boundaries,  and  obtain  also  for  the  case 
of  intrinsic  distance  functions  on  sub-manifolds  of  ]Rd ,  a  computationally  optimal  approach.  For  point 
clouds,  the  offset  band  is  constructed  without  the  need  to  explicitly  find  the  underlying  manifold,  thereby 
computing  intrinsic  distance  functions  and  geodesics  on  point  clouds  while  skipping  the  manifold  recon¬ 
struction  step.  The  case  of  point  clouds  representing  noisy  samples  of  a  sub-manifold  of  Euclidean 
space  is  studied  as  well.  All  the  underlying  theoretical  results  are  presented,  together  with  experimental 
examples,  and  comparisons  to  graph-based  distance  algorithms. 


1  Introduction 

One  of  the  most  popular  sources  of  point  clouds  are  3D  shape  acquisition  devices,  such  as  laser  range 
scanners,  with  applications  in  geoscience,  art  (e.g.,  archival),  medicine  (e.g.,  prohestetics),  manufacturing 
(from  cars  to  clothes),  and  security  (e.g.,  recognition),  among  other  disciplines.  These  scanners  provide  in 
general  raw  data  in  the  form  of  unorganized  point  clouds  representing  surface  samples.  With  the  increasing 
popularity  and  very  broad  applications  of  this  source  of  data,  it  is  natural  and  important  to  work  directly 
with  this  representation,  without  having  to  go  to  the  intermediate  step  of  fitting  a  surface  to  it  (step  that  can 
add  computational  complexity  and  introduce  errors).  See  for  example  [6,  9,  10,  12,  13,  20,  21,  24,  25]  for 
a  few  of  the  recent  works  with  this  type  of  data.  Note  that  point  clouds  can  also  be  used  as  primitives  for 
visualization,  e.g.,  [7,  13,  26],  as  well  as  for  editing  [34]. 

Another  important  field  where  point  clouds  are  found  is  in  the  representation  of  high-dimensional  man¬ 
ifolds  by  examples  (see  for  example  [15,  19,  30]).  This  type  of  high-dimensional  and  general  co-dimension 
data  appeal's  in  almost  all  disciplines,  from  computational  biology  to  image  analysis  to  financial  data. 

Note  that  in  general  a  point  cloud  representation  is  dimension  free,  in  contrast  with  other  popular  rep¬ 
resentations  such  as  triangular  meshes.  Some  operations,  such  as  the  union  of  point  clouds  acquired  from 
multiple  views,  are  much  easier  than  when  performed  on  the  triangular  meshes  obtained  from  them.  This 
paper  addresses  one  of  the  most  fundamental  operations  in  the  study  and  processing  of  sub-manifolds  of 
Euclidean  space,  the  computation  of  intrinsic  distance  functions  and  geodesics.  We  show  that  this  can  be 
done  directly  working  with  the  point  cloud,  without  the  need  for  reconstructing  the  underlying  manifold.  We 
present  the  corresponding  theoretical  results,  experimental  examples,  and  basic  comparisons  to  mesh-based 
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distance  algorithms.  The  results  arc  valid  for  general  dimensions  and  co-dimensions,  and  for  manifolds  with 
or  without  boundary.  These  results  include  the  analysis  of  noisy  point  clouds  obtained  from  sampling  the 
manifold.  The  theoretical  results  arc  presented  without  proof  due  to  space  limitations.  The  corresponding 
proofs  can  be  found  in  our  extended  report  [23]. 

A  number  of  key  building  blocks  form  paid  of  the  framework  here  introduced.  The  first  one  is  based  on 
the  fact  that  distance  functions  intrinsic  to  a  given  sub-manifold  of  Fid  can  be  accurately  approximated  by 
Euclidean  distance  functions  computed  in  a  thin  offset  band  that  surrounds  this  manifold.  This  concept  was 
first  introduced  in  [22],  where  convergence  results  arc  given  for  co-dimension  one  hyper-surfaces  without 
boundary.  This  result  is  reviewed  in  §2.  In  this  paper,  we  first  extend  this  to  general  co-dimension  and  to 
deal  with  manifolds  with  or  without  boundary,  §3.  We  also  show  that  the  approximation  is  true  not  only 
for  the  intrinsic  distance  function  but  also  for  the  intrinsic  geodesic.  This  is  not  a  straightforward  corollary, 
since  geodesics  arc  based  on  the  gradient  of  the  distance  function, 1  which  contains  singularities  at  the  cut 
locus  [33].  The  approximation  of  intrinsic  distance  functions  (and  geodesics)  by  extrinsic  Euclidean  ones 
permits  to  compute  them  using  computationally  optimal  algorithms  in  Cartesian  grids.  These  algorithms 
arc  based  on  the  fact  that  the  distance  function  satisfies  a  Hamilton-Jacobi  partial  differential  equation  (see 
§2),  for  which  consistent  and  fast  algorithms  have  been  developed  in  Cartesian  grids  [14,  28,  29,  32] 2  (see 
[16]  for  extensions  to  triangular  meshes  and  [31]  for  other  Hamilton-Jacobi  equations).  That  is,  due  to  these 
results,  we  can  use  computationally  optimal  algorithms  in  Cartesian  grids  (with  boundaries)  also  to  compute 
distance  functions,  and  from  them  geodesics,  intrinsic  to  a  given  manifold,  and  in  a  computationally  optimal 
fashion.  Note  that  in  contrast  with  the  popular-  Dijkstra  algorithm,  these  numerical  techniques  are  consistent, 
they  converge  to  the  true  distance  when  the  grid  is  refined. 

Once  this  basic  results  are  available,  we  can  then  proceed  and  work  with  point  clouds.  The  basic  idea 
here  is  to  construct  the  offset  band  directly  from  the  point  cloud  and  without  the  intermediate  step  of  man¬ 
ifold  reconstruction.  This  is  addressed  in  §4  and  §5  for  noise-free  points  and  manifold  samples,  and  in  §6 
for  points  considered  to  be  noisy  samples  of  the  manifold.  In  this  case,  we  explicitly  compute  the  proba¬ 
bility  that  the  constructed  offset  band  contains  the  underlying  manifold.  As  we  expect,  this  probability  is 
a  function  of  the  number  of  point  samples,  the  noise  level,  the  size  of  the  offset,  and  the  basic  geometric 
characteristics  of  the  underlying  manifold.  In  the  experimental  section,  §7,  we  also  show  the  advantages  of 
this  framework  over  mesh-based  distance  computations. 

To  conclude  this  introduction,  we  should  note  that  to  the  best  of  our  knowledge,  the  only  additional 
work  explicitly  addressing  the  computation  of  distance  functions  and  geodesics  for  point  clouds  is  the  one 
reported  in  [4], 3  The  comparison  of  our  framework  with  the  one  proposed  by  this  paper  is  deferred  to  §8, 
where  we  also  report  the  directions  our  research  is  taking. 

2  Preliminary  results  and  notation 

In  this  section  we  briefly  review  the  main  results  in  [22],  where  the  idea  of  approximating  intrinsic  distances 
and  geodesics  by  extrinsic  ones  was  first  introduced. 

We  first  introduce  some  basic  notation  that  will  be  used  throughout  the  article.  For  a  compact  and 
connected  set  A  E  Md,  cIa( •,•)  denotes  the  intrinsic  distance  between  any  two  points  of  A,  measured 

’Geodesics  are  the  integral  curves  corresponding  to  the  gradient  directions  of  the  intrinsic  distance  function,  and  are  obtained 
back-propagating  in  this  gradient  direction  from  the  target  point  to  the  source  one. 

“Tsitsiklis  first  described  an  optimal-control  type  of  approach  to  solve  the  Hamilton-Jacobi  equation,  while  independently 
Sethian  and  Helmsen  both  developed  techniques  based  on  upwind  numerical  schemes. 

’In  addition  to  studying  the  computation  of  distance  functions  on  point  clouds,  [4,  30]  address  the  important  combination  of  this 
with  multidimensional  scaling  for  manifold  analysis.  Prior  work  on  using  geodesics  and  multidimensional  scaling  can  be  found  in 
[27], 


2 


by  paths  constrained  to  be  in  A.  Given  a  A; -dimensional  sub-manifold  M.  of  Md,  denotes  the  set 

{. x  E  lRd  :  d(A4,x)  <  h}.  This  is  basically  an  h- offset  of  Ad.  To  state  that  the  sequence  of  functions 

n 

{/n(')}npi+  uniformly  converges  to  /(•)  as  n  t  oo,  we  frequently  write  fn  =4  /.  For  a  given  event  £, 
P  (£)  stands  for  its  probability  of  occurring.  For  a  random  variable  (R.V.  from  now  on)  X,  its  mean  value 
is  denoted  by  E  (X).  For  a  function  f  :  £1  M,  and  a  subset  A  of  Q,  /  |  4  :  A  — >  FI  denotes  the  restriction 
of  /  to  A. 

In  [22],  we  presented  a  new  approach  for  the  computation  of  weighted  intrinsic  distance  functions  over 
hyper-surfaces.  We  proved  convergence  theorems  and  also  addressed  the  fast,  computationally  optimal, 
computation  of  such  approximations.  The  key  starting  idea  is  that  distance  functions  satisfy  the  (intrinsic) 
Eikonal  equation,  a  particular  case  of  the  general  class  of  Hamilton- Jacobi  partial  differential  equations. 
Given  p  e  S  (an  hyper-surface  in  Md),  we  want  to  compute  ds(p,  ■)  :  S  — >  1R+  U  {0},  the  intrinsic  distance 
function  from  every  point  on  S  to  p.  It  is  well  known  that  the  distance  function  ds(p,  •)  satisfies,  in  the 
viscosity  sense,  the  equation 

f  ||V5<is(p,s)||  =  1  V*e<S 

1  ds(p,p)  =  0 

where  V5  is  the  intrinsic  differentiation  (gradient).  Instead  of  solving  this  intrinsic  Eikonal  equation  on  S , 
we  solve  the  corresponding  extrinsic  one  in  the  offset  band  £1$: 

f  l|V*dnh(p,x)||  =  1  Vx  E 

\  dnh(p,p)  =  0 

where  d<}h  (p.  •)  is  the  Euclidean  distance  and  therefore  now  the  differentiation  is  the  usual  one.  The  main 
result  of  [22]  is 


Theorem  1  ([22])  Let  p  and  q  be  any  two  points  on  the  smooth  hyper-surface  S,  then 
ds  (p,  <l)  ~  d^h  (p,  q)  I  <  Csfh,  for  small  enough  h,4  where  C  is  a  constant. 


This  simplification  of  the  intrinsic  problem  into  an  extrinsic  one  permits  the  use  of  the  computationally 
optimal  algorithms  mentioned  in  the  introduction.  This  makes  computing  intrinsic  distances,  and  from  them 
geodesics,  as  simple  and  computationally  efficient  as  computing  them  in  Euclidean  spaces.  Moreover,  as 
detailed  in  [22],  the  approximation  of  the  intrinsic  distance  ds  by  the  extrinsic/Euclidean  one  dnh  is  never 
less  accurate  than  the  numerical  error  of  these  algorithms. 

In  [22],  the  result  above  was  limited  to  co-dimension  one  closed  hyper-surfaces,  and  the  theory  was 
applied  to  implicit  surfaces,  where  computing  the  offset  band  is  straightforward.  It  is  the  purpose  of  the 
present  work  to  extend  the  aforementioned  Theorem  to  deal  with:  (1)  sub-manifolds  of  lRd  of  any  codi¬ 
mension  and  with  boundary,  (2)  sub-manifolds  of  Md  represented  as  points  cloud,  (3)  random  sampling  of 
sub-manifolds  of  Md  in  presence  or  absence  of  noise,  and  (4)  convergence  of  geodesic  curves  in  addition  to 
distance  functions.  We  should  note  that  Theorem  1  holds  even  when  the  metric  is  not  the  one  inherited  from 
Md ,  obtaining  weighted  distance  functions,  see  [22].  Although  we  will  not  present  these  new  results  in  such 
generality,  this  is  a  simple  extension  that  will  be  reported  elsewhere. 


3  General  distance  functions  and  geodesics  approximation  results 

We  now  generalize  Theorem  1  to  cover  open  manifolds  and  general  co-dimension.  The  convergence  of  the 
corresponding  geodesic  curves,  which  was  not  studied  in  [22],  is  obtained  as  well.  As  mentioned  above, 

4“Small  enough  h”  means  that  h  <  1/  max*  Ki(S).  where  k,(S)  is  the  i-th  principal  curvature  of  S.  This  guarantees  having 
smoothness  in  see  [22], 
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proofs  arc  omitted  due  to  space  constraints,  and  they  can  be  found  in  [23]. 

Theorem  2  Let  S  be  a  compact  C  ’1  sub-manifold  of  Md  with  smooth  boundary  dS.  We  then  have 

/40 

1.  Uniform  convergence  of  the  distance  functions:  d^h  |5x5(-,-) 

2.  Convergence  of  the  geodesics:  Let  x  and  y  <E  S  he  joined  by  a  unique  minimizing  geodesic  7  s  ■  [0,1]  ^<5 

/40 

over  S,  and  let  y^  :  [0, 1]  -»  be  the  corresponding  minimizing  geodesic  in  the  band  Llg.  Then  y^  =t  7 5. 

Note  that  the  new  setting  is  wider  than  the  one  in  Theorem  1  from  [22],  since  the  co-dimension  of  the 
underlying  manifold  is  not  necessarily  1.  This  is  very  important  for  applications  such  as  dimensionality 
reduction,  where  the  dimension  of  the  underlying  manifold  is  unknown  beforehand.  The  convergence  of 
the  distances  is  uniform  following  the  same  hypotheses  as  in  Theorem  1,  but  we  have  forfeited  the  exact 
rate  of  convergence  estimates.  To  obtain  these  rates  for  open  manifolds  of  general  co-dimension,  further 
restrictions  arc  imposed  in  dS,  see  Theorem  3. 

We  can  immediately  obtain  the  following  corollary  which  will  become  useful  when  we  need  to  shape 
the  Ll[ j  band  in  more  convenient  ways. 

Corollary  1  Let  S  and  dS  satisfy  the  hypotheses  of  Theorem  2.  Let  { be  a  family  of  sets  in 

Md  such  that  S  C  E,  V  %  €  IN  and  d%(Ej,<S)  q  stands  for  the  Hausdorff  distance).  Then 

tt+oo 

|sxs  (-,  ■)  =3 

We  now  restrict  the  manifold  in  order  to  obtain  explicit  rates  of  convergence  for  general  co-dimension 
open  manifolds. 

Definition  1  We  say  that  the  manifold  S  with  boundary  dS  is  strictly  convex  if  for  every  pair  of  points  x 

O  O 

and  y  in  the  interior  S  of  S,  there  exists  a  unique  geodesic  joining  them  that  is  entirely  contained  in  S- 


Theorem  3  (Rate  of  Convergence)  Under  the  hypotheses  of  Theorem  2,  and  assuming  dS  to  be  either 
void  or  strictly  convex,  we  have  that  for  small  enough  h  >  0  and  a  constant  C  that  does  not  depend  on  h, 


dnh(x,y) 


ds{x,y )  <  Cy/hds{x,y),  and  then  ma <tx.(x,y)eSxS 


dnh{x,y) 


ds{x,y ) 


<  CVh. 


Let’s  point  out  that  we  do  not  know  yet  whether  the  convexity  condition  for  dS  is  really  necessary  in 
order  to  prove  this  rate  of  convergence,  and  advances  on  this  precise  issue  will  be  reported  elsewhere.  Note 
that  still,  from  Theorem  2,  we  have  uniform  convergence,  even  if  the  rate  is  not  explicitly  computed. 

To  summarize  this  section,  the  extrinsic  Euclidean  distance  and  geodesic  computed  in  a  thin  offset  band 
around  the  manifold  S  uniformly  converges  to  the  intrinsic  ones  on  <S,  for  general  co-dimensions  and  both 
closed  and  open  manifolds.  Under  special  cases,  the  rate  of  convergence  can  be  explicitly  computed.  With 
these  new  results,  we  arc  now  ready  to  proceed  with  manifolds  given  by  point  clouds. 


4  Point  clouds 

Let  Vn  =  \pi , . . . , pTl }  be  points  on  the  sub-manifold  S  and  define  Qtp  =  (JILi  B(Pi>h).  Let  also  h 
and  Vn  be  such  that  S  C  0^  ,5  Based  on  this,  we  then  have  S  C  C  f Vf.  Now  the  argument 

carries  over  easily,  for  if  p  and  q  arc  any  two  points  in  S,  one  has  dnh  (p,  q)  <  dnh  (p.  q)  <  ds(p,  9),  so 

<5  Vn 

0  <  ds  (p,  q)  —  d(Ah  [p,  q )  <  ds  (p,  q)  —  dnh  (p,  q),  and  the  rightmost  quantity  can  be  bounded  by  C  h 1  ;/'2 

Vn 

(see  [22]  and  §3).  We  therefore  obtain  the  following  basic  results  for  point  clouds: 

3We  should  point  out  that  results  such  as  those  in  [1,  2,  3]  provide  basic  theory  to  guarantee  that  the  offset  band  ,  while 
covering  S,  has  the  correct  topology. 
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Theorem  4  Assume  that  S  C  ilf  then,  under  the  hypotheses  of  Theorem  3,  we  have  that  for  a  constant  C, 


max, 


(p,q)eSxS 


ds{p,q)  -  dQh  (p,  q) 

rn 


<  CVh. 


Corollary  2  Let  { hn  }neiN  be  such  that  hn  - — >  0  and  S  C  {Vfl  for  all  n,  then  dQhn  \sxS  converges 

n  S  LVn 

uniformly  to  ds  asnf  +oo. 

These  basic  results  indicate  that  in  order  to  compute  distance  functions,  and  from  them  geodesics,  on 
manifolds  represented  by  point  clouds,  all  what  needs  to  be  done  is  to  construct  the  union  of  balls  f l  If,  ,  and 
then  inside  this  band  in  lRd  use  classical  Euclidean  algorithms.  Note  than  the  explicit  reconstruction  of  the 
underlying  manifold  is  not  needed. 

In  general,  different  radii  for  each  ball  defining  the  set  f lip  can  be  considered.  In  this  case  let  VnlZn  be 
given  by  {(pi,  ri), . . .  ,  (pn ,  rn ) } .  that  is  the  set  of  pairs  (p.  r )  where  p  arc  points  sampled  on  the  manifold. 

Define  =  (J"=1  B(pi,  rf).  Then,  if  we  ensure  that  S  C  Op"  and  that  R  =  max i <.,;<«  N  is  small 
enough,  we  can  prove  results  similar  to  Theorem  4  and  Corollary  2.  In  general,  all  what  is  needed  is 
to  guarantee  that  the  union  of  balls  includes  the  underlying  manifold  S  (this  coverage  principle  will  be 
addressed  in  more  detail  in  future  sections).  Moreover,  we  arc  not  restricted  to  use  just  balls,  and  for 
example  hyper-ellipsoids  can  be  used  as  well  to  define  the  band.  In  particular,  we  can  use  hyper-ellipsoids 
that  arc  adapted  to  the  local  geometry  and  density  of  the  point  cloud.  This  can  be  obtained  by  iterative  local 
principal  component  analysis  [8],  This  scenario  may  be  desirable  for  example  when  more  resources  arc 
spent  on  areas  of  the  manifold  that  arc  “harder”  (high  curvatures)  than  on  those  geometrically  simpler  parts 
(low  curvatures). 

5  Distance  functions  on  randomly  sampled  manifolds 

In  the  previous  section,  we  dealt  with  a  deterministic  set  of  points  in  Ft'1.  We  now  consider  a  probabilistic 
framework,  with  points  randomly  sampled  from  a  manifold  S  in  IRd.  As  one  can  guess  from  §4,  we  first 
have  to  study  the  crucial  event  {5  C  0^},  when  we  arc  given  that  the  points  in  Vn  arc  randomly  sampled 
from  S.  To  illustrate  our  results,  we  assume  a  uniform  distribution  of  the  points  over  the  manifold,  Vn  is  an 
i.i.d.  sequence  of  points  uniformly  sampled  over  S. 6  This  means  that  for  any  A  C  S,P  (pi  E  A)  =  fygj, 
for  all*  E  {1,2,...  .  n}.  Here,  p  (A)  denotes  the  volume  measure  of  the  subset  A  of  the  manifold  S.  The 
results  below  can  be  extended  for  other  sampling  strategies.7 

Theorem  5  Let  S  be  a  smooth  compact  k-dimensional  sub-manifold  of  LRd.  Let  n  t  +oo  and  hn  j,  0  be 
such  that  both  limn  =  0  and  limn  nh *  =  +oc  hold,  then  p|<SC  0^,"  j  ^ 

Of  independent  interest  is  the  following  Random  Sampling  Asymptotic  Theorem  for  sub-manifolds  of 
Fl'K  which  can  be  proven  using  the  preceding  Theorem:8 

Theorem  6  Under  the  hypotheses  of  Theorem  5,  we  have  plimn  dy  (S.  )  =  0. 

We  now  proceed  with  the  convergence  results  for  the  corresponding  distance  functions.  Being  p  and  q 
points  on  S,  our  main  goal  is  to  prove  some  kind  of  random  convergence  towards  ds(p,  q)  for  the  random 

6 We  will  also  write  p,  ~  U  [<!>]. 

7Note  that  we  are  considering  a  sequence  of  offsets  that  goes  to  0.  This  is  not  necessary,  but  it  is  in  the  spirit  of  our  work,  since 
we  need  hn  — >  0  in  order  to  achieve  convergence  of  the  distances. 

8Here,  for  the  sequence  of  R.V.  Xn,  and  the  R.V.  X,  p  limn  Xn  =  X  means  that  for  any  5  >  0,  limn  P  ( \Xn  —  X\  >  <5)  =  0. 
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variables  d(lh  (p,  q)  .  Let’s  assume  from  now  on  that  h  (or  hn )  is  small  enough.  In  the  same  spirit  as  in 

'Pn 

the  deterministic  case,  it  is  meaningful  to  consider  max^  \dQh  (x.  y )  —  ds(x.  y)\.  However,  to  this 

end,  we  must  be  able  to  evaluate  the  function  dnh  (•,  ■)  at  all  points  on  S  x  S.  This  can  be  guaranteed  only 

Pn 

if  S  C  flip  .  We  are  then  led  to  consider  the  event  in  the  following  Theorem: 

Theorem  7  Under  the  hypotheses  of  Theorem  5,  and  for  e  >  0, 

lim^ooP  ({{<S  CSlJrJn  {max^esxs  \dQhn  (x,y)  -  ds(x,y)\  <  £>})  =  1- 

Details  about  how  the  intrinsic  characteristics  of  the  manifold  play  a  role  in  the  above  convergence 
results  arc  reported  in  [23], 

To  summarize,  in  this  section  we  have  addressed  the  distance  convergence  results  for  randomly  sampled 
manifolds.  We  have  presented  results  regarding  the  probability  that  the  manifold  of  interest  is  covered  by 
the  union  of  balls  centered  at  the  random  samples  and  that  the  error  between  the  corresponding  intrinsic  and 
extrinsic  distance  functions  is  bounded.  We  now  proceed  to  add  noise  to  the  samples. 

6  Distance  functions  on  noisily  sampled  manifolds 

Since  manifold  samples  arc  expected  to  be  noisy,  it  is  important  to  add  this  to  the  results  presented  above. 
Although  a  particular  model  of  noise  is  assumed  in  this  section,  the  results  can  be  easily  extended  to  other 
models. 

Assume  that  {pi, . . .  ,pn}  is  a  set  of  i.i.d.  random  points  such  that  each  pt  ~  U[flA].  The  main  point 
here  is  to  notice  that,  whenever  {5  C  }  holds,  we  have: 

max  | dnh  {x,y)  -  ds(x,y)\  <  Cs/h  +  A.  (1) 

(x,y)eSxS  vn 

In  fact,  as  a  consequence  of  the  presence  of  noise  in  the  points,  we  have  that  ma xpe-pn  d(S,p)  <h  +  A. 
Therefore,  Llfpn  C  Ll^+A,  and  the  claim  follows  from  Theorem  4.  That  is,  our  bound  regarding  the  absolute 
convergence  of  the  distance  in  the  cloud  of  points  is  worse  than  the  one  we  can  guarantee  in  absence  of 
noise.  Moreover,  we  have  bounded  this  “worsening”  using  the  parameter  describing  the  noise.  Up  to  first 
order,  this  worsening  can  be  approximated  by  A4=. 

To  conclude,  we  wish  to  estimate  the  probability  of  having  S  C  (J'l_  (  B  (pt .  h ).  It  is  easy  to  see  that  as  a 
first  “reality  compliant”  condition  one  should  have  h  >  A.  Starting  from  the  setting  for  Theorem  7,  we  can 
prove  the  following: 

Theorem  8  Let  S  be  a  k-dimensional  sub-manifold  of  lRd  satisfying  the  usual  hypotheses.  Let  An  < 
hn  — >  0  as  n  t  +oo,  such  that  both  nAj  — >  oo  and  ^  -»  0  as  n  t  +oo  hold,  then  for  any  e  >  0 

lim^ooP  ({{<S  C  Up”  }  n  {max(a.;3/)e<Sx<5  \dQhn  (x,y)  -  ds(x,y)\  <  e}})  =  1- 

7  Examples 

We  now  present  examples  of  distance  functions  and  geodesics  for  point  clouds.  Figure  1;  use  these  compu¬ 
tations  to  find  intrinsic  Voronoi  diagrams.  Figure  2  (see  also  [17,  18]);  and  compare  the  results  with  those 
obtained  with  mesh-based  techniques.  Figure  3. 9  The  theoretical  results  presented  in  previous  sections  show 

9 All  the  figures  in  this  paper  are  in  color.  VRML  files  corresponding  to  these  examples  can  be  found  at 
mountains.ece.umn.edu/~guille/pc.htm. 
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that  the  intrinsic  distance  and  geodesics  can  be  approximated  by  the  Euclidean  ones  computed  in  the  band 
defined  by  the  union  of  balls  centered  at  the  points  of  the  cloud.  The  problem  is  then  simplified  to  first  com¬ 
puting  this  band,  and  then  use  computationally  optimal  techniques  to  compute  the  distances  and  geodesics 
inside  this  band,  exactly  as  done  in  [22]  for  implicit  surfaces.  The  band  itself  can  be  computed  in  several 
ways,  and  for  the  examples  below  we  have  used  constant  radii.  Locally  adaptive  radii  can  be  obtained  as 
mentioned  in  §4,  or  can  be  based  on  diameters  obtained  for  example  from  minimal  spanning  trees. 


Figure  1:  Left:  Intrinsic  distance  function  for  a  point  cloud.  A  point  is  selected  in  the  head  of  the  David,  and  the 
intrinsic  distance  is  computed  following  the  framework  here  introduced.  The  point  cloud  is  colored  according  to  their 
intrinsic  distance  to  the  selected  point,  going  from  dark  blue  (close)  to  bright  red  (far).  The  offset  band,  given  by  the 
union  of  balls,  is  shown  next  to  the  distance  figure.  Right:  Same  as  before,  with  a  geodesic  curve  between  two  selected 
points. 


Figure  2:  Voronoi  diagram  for  point  clouds.  Four  points  (left)  and  two  points  (right)  are  selected  on  the  cloud,  and 
the  point  cloud  is  divided  ( colored)  according  to  their  geodesic  distance  to  these  four  points.  Note  that  this  is  a  surface 
Voronoi,  based  on  geodesics  computed  with  our  proposed  framework,  not  an  Euclidean  one. 

We  now  make  some  very  basic  comparisons  between  our  approach  to  geodesic  distance  computations 
and  those  based  on  graph  approximations  to  the  manifold,  such  as  the  one  in  Isomap  [30]. 10  The  goal  is 
to  show  that  such  graph-based  techniques  are  more  sensitive  to  noise  in  the  point  cloud  sample.  This  is 
expected,  since  the  geodesic  in  such  techniques  goes  through  the  noisy  samples,  while  in  our  approach,  they 
just  go  through  the  union  of  balls. 

In  the  case  of  N  noisy  points  (where  the  noise  is  as  in  the  hypotheses  of  Section  §6,  with  A  -C  1), 
sampled  on  a  rectilinear  segment  of  length  do  in  M2,  an  explicit  construction  shows  that  the  mean  value  error 
for  the  graph  approximation  df)  is  approximately  given  by  E  {\df";  —  do  | )  —  On  the  other  hand,  for 

our  approximation,  dft,  if  the  segment  is  contained  in  the  union  of  the  balls  centered  at  the  sampling  points, 
10Isomap  builds  a  mesh  by  locally  connecting  the  (noisy)  samples. 
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E  [\d^  —  do|)  =0.  The  probability  of  covering  the  segment  by  the  band  can  be  made  arbitrarily  close  to 
1  by  increasing  N.  More  precisely,  one  can  prove  that  if  p  stands  for  the  value  of  the  probability  of  not 
covering  the  segment,  then  p  <  k^(  1  —  k'-^)N,  for  some  positive  constants  k  and  k' .  We  then  note  that  for 
a  fixed  noise  level  A,  by  increasing  N  we  actually  worsen  the  graph  approximation,  whereas  we  are  making 
our  approximation  better. 

In  the  table  below  we  present  results  of  simulations  earned  out  for  the  SwissRoll  dataset  [30],  see  Fig¬ 
ure  3  .  We  used  10000  points  to  define  the  manifold.  We  then  generated  10000  noise  vectors,  being  each 
component  uniform  with  power  one  and  zero  mean.  Then,  we  generated  noisy  datasets  from  the  noise¬ 
less  SwissRoll  dataset  by  adding  the  noise  vector  times  a  constant  n *  to  each  vector  of  the  noiseless  initial 
dataset.  We  then  chose  1000  corresponding  points  in  each  dataset  and  computed  the  intrinsic  pairwise  dis¬ 
tance  approximation  obtaining  the  matrices  {(.D??"*)}  and  { (Dl-jrik ) }  for  the  graph-based  and  our  approach 
respectively,  where  k  =  1, 2, . . .  ,  5,  i,  j  6  [1, 1000],  and  denotes  the  noise  level.  We  then  computed  the 
values  of  rnax.lJ\D^nk  —  Dfj0 |  and  max.l3  \D>-':'lk  —  D^j°\  for  each  k,  where  D!k(>  and  D^j°  stand  for  noise¬ 
less  intrinsic  distance  approximations.  In  the  table,  h  indicates  the  radii  and  k  the  size  of  the  neighborhood 
for  Isomap.  The  graph  approximation  shows  less  robustness  to  noise  than  our  method,  as  was  argued  above. 
This  is  also  true  for  the  sensitivity, 1 1  where  our  approach  outperforms  the  graph-based  one  by  at  least  one 
order  of  magnitude.  Note  that  the  sensitivity  for  our  approach  can  be  formally  studied  from  Theorem  3. 


Noise  Power  (n\) 

maxi3\D\r  ~  mf\ 

k 

maxifD^r  ~  D^\ 

h 

0.0001 

2.5222 

7 

0.5266 

1.8 

0.01 

4.6409 

7 

0.9430 

1.8 

0.04 

5.1737 

7 

1.2489 

1.8 

0.09 

5.3292 

7 

1.4682 

1.8 

0.16 

5.4651 

7 

1.7965 

1.8 

Figure  3:  Surface  used  to  compare  graph-based  approaches  with  our  framework. 


8  Concluding  remarks 


In  this  paper,  we  have  shown  how  to  compute  distance  functions  and  geodesics  intrinsic  to  a  generic  manifold 
defined  by  a  point  cloud,  without  the  intermediate  step  of  manifold  reconstruction.  The  basic  idea  is  to 
use  computationally  optimal  algorithms  for  computing  Euclidean  distances  in  an  offset  band  surrounding 
the  manifold,  and  use  these  to  approximate  the  intrinsic  distance.  The  underlying  theoretical  results  were 
complemented  by  experimental  illustrations. 

As  mentioned  in  the  introduction,  an  alternative  technique  to  compute  geodesic  distances  was  introduced 
in  [4,  30].  In  contrast  with  our  work,  the  effects  of  noise  were  not  addressed  in  [4],  Moreover,  as  we  have 
explicitly  shown,  our  framework  is  more  robust  to  noise.  We  should  also  add  that  the  theoretical  results 
reported  in  [4]  impose  stronger  restrictions  on  the  manifolds  than  the  ones  we  need.  On  the  other  hand,  we 


i  „  ...  .  ,  -  ,  L  distance  for  noisy  points 

Sensitivity  is  denned  as  1  —  -t-— - T - J  r  ■  . 

J  distance  for  clean  points 


should  note  that  the  memory  requirements  of  our  framework  arc  larger,  and  this  needs  to  be  addressed  for 
very  high  dimensions  [23]. 

We  arc  currently  working  on  the  use  of  this  framework  to  create  multiresolution  representations  of  point 
clouds  (in  collaboration  with  N.  Dyn,  see  also  [6,  9,  10,  25]),  to  perform  object  recognition  (see  [11]  for 
early  work  on  the  use  of  geodesic  distances  for  this),  and  to  compute  basic  geometric  characteristics  of  the 
underlying  manifold,  all  this  of  course  without  reconstructing  the  manifold.  Applications  of  our  framework 
for  high  dimensional  data  arc  also  currently  being  addressed,  and  preliminary  results  arc  reported  in  [23]. 


Acknowledgments 

We  acknowledge  useful  conversations  on  the  topic  of  this  paper  with  L.  Aspirot,  P.  Bermolen,  R.  Coifman, 
N.  Dyn,  O.  Gil,  R.  Kimmel,  and  A.  Pardo.  We  thank  M.  Levoy  and  the  Digital  Michelangelo  Project  for 
data  provided  for  this  project.  This  work  was  supported  by  a  grant  from  the  Office  of  Naval  Research 
ONR-NOOO 14-97- 1-0509,  the  Presidential  Early  Career  Award  for  Scientists  and  Engineers  (PECASE),  and 
a  National  Science  Foundation  CAREER  Award. 


References 

[1]  N.  Amenta,  S.  Choi,  and  R.  Kolluri,  “The  power  crust,  unions  of  balls,  and  the  medial  axis  transform,”  Computa¬ 
tional  Geometry:  Theory  and  Applications  19,  pp.  127-153,  2001. 

[2]  N.  Amenta,  S.  Choi,  and  R.  Kolluri,  “The  power  crust,”  Proceedings  of  6th  ACM  Symposium  on  Solid  Modeling, 
pp.  249-260,  2001. 

[3]  N.  Amenta  and  R.  Kolluri,  “Accurate  and  efficient  unions  of  balls,”  ACM  Symposium  on  Computational  Geometry, 
pp.  119-128,2000. 

[4]  M.  Bernstein,  V.  de  Silva,  J.  Langford,  and  J.  Tennenbaum,  “Graph  approximations  to  geodesics  on  embeddeed 
manifolds,”  http  :  //isomap  .  Stanford.  edu/BdSLT  ,ps 

[5]  Blitz++  website:  www .  oonumerics  .  org/blitz 

[6]  J-D.  Boissonnat  and  F.  Cazals,  “Coarse-to-fine  surface  simplification  with  geometric  guarantees,”  in  A.  Chalmers 
and  T-M.  Rhyne,  Editors,  EUROGRAPHICS  ’01,  Manchester,  2001. 

[7]  M.  Botsch,  A.  Wiratanaya,  and  L.  Kobbelt,  “Efficient  high  quality  rendering  of  point  sampled  geometry,”  EURO¬ 
GRAPHICS  Workshop  on  Rendering,  2002. 

[8]  R.  Coifman,  personal  communication. 

[9]  T.  K.  Dey,  J.  Giesen,  and  J.  Hudson,  “Decimating  samples  for  mesh  simplification,”  Proc.  13th  Canadian  Confer¬ 
ence  on  Computational  Geometry,  pp.  85-88,  2001. 

[10]  N.  Dyn,  M.  S.  Floater,  and  A.  Iske,  “Adaptive  thinning  for  bivariate  scattered  data,”  Journal  of  Computational 
and  Applied  Mathematics  145(2),  pp.  505-517,  2002. 

[11]  A.  Elad  (Elbaz)  and  R.  Kimmel,  “Bending  invariant  representations  for  surfaces,”  In  Proc.  ofCVPR’01,  Hawaii, 
Dec.  2001. 

[12]  M.  S.  Floater  and  A.  Iske,  “Thinning  algorithms  for  scattered  data  interpolation,”  BIT  Numerical  Mathematics 
38,  pp.  705-720, 1998. 

[13]  M.  Gross  et  al.  Point  Based  Computer  Graphics,  EUROGRAPHICS  Lecture  Notes,  2002 
(graphics.stanford.edu/~niloy/research/papers/ETH/PointBasedComputerGraphics_TutorialNotes.pdf). 


9 


[14]  J.  Helmsen,  E.  G.  Puckett,  P.  Collela,  and  M.  Dorr,  “Two  new  methods  for  simulating  photolithography  develop¬ 
ment  in  3D,”  I1  roc.  SPIE  Microlithography  IX,  pp.  253,  1996. 

[15]  P.  W.  Jones,  “Rectifable  sets  and  the  traveling  salesman  problem,”  Invent.  Math.  102,  pp.  1-15,  1990. 

[16]  R.  Kimmel  and  J.  A.  Sethian,  “Computing  geodesic  paths  on  manifolds,”  Proc.  National  Academy  of  Sciences 
95:15,  pp.  8431-8435, 1998. 

[17]  R.  Kunze,  F.  E.  Wolter,  and  T.  Rausch,  “Geodesic  Voronoi  diagrams  on  parametric  surfaces,”  IEEE ,  1997. 

[18]  G.  Leibon  and  D.  Letscher,  “Delaunay  triangulations  and  Voronoi  diagrams  for  Riemannian  manifolds,”  Com¬ 
putational  Geometry  2000 ,  Hong  Kong,  2000. 

[19]  G.  Lerman,  “How  to  partition  a  low-dimensional  data  set  into  disjoint  clusters  of  different  geometric  structure,” 
preprint,  2000. 

[20]  L.  Linsen,  “Point  cloud  representation,”  CS  Technical  Report,  University  of  Karlsruhe,  2001. 

[21]  L.  Linsen  and  H.  Prautzsch,  “Local  versus  global  triangulations,”  EUROGRAPHICS  ’01,  2001. 

[22]  F.  Memoli  and  G.  Sapiro,  “Fast  computation  of  weighted  distance  functions  and  geodesics  on  implicit  hyper¬ 
surfaces,”  Journal  of  Computational  Physics  173:2,  pp.  730-764, 2001. 

[23]  F.  Memoli  and  G.  Sapiro,  “Distance  functions  and  geodesics  on  surfaces  and  point  clouds,”  January  2003 
(mountains.ece.umn.edu/~guille/pc.htm). 

[24]  M.  Pauly  and  M.  Gross,  “Spectral  processing  of  point-sampled  geometry,”  ACM  SIGGRAPH,  pp.  379-386, 2001 

[25]  M.  Pauly,  M.  Gross,  and  L.  Kobbelt,  “Efficient  simplification  of  point-sampled  surfaces,”  IEEE  Visualization 

2002. 

[26]  S.  Rusinkiewicz  and  M.  Levoy,  “QSplat:  A  multiresolution  point  rendering  system  for  large  meshes,”  Computer 
Graphics  (SIGGRAPH  2000  Proceedings),  2000. 

[27]  E.  Schwartz,  A.  Shaw,  and  E.  Wolfson,  “A  numerical  solution  to  the  generalized  mapmaker’s  problem:  Flattening 
nonconvex  polyhedral  surfaces,”  IEEE  Transactions  on  Pattern  Analysis  and  Machine  Intelligence  11:9,  pp.  1005 
-1008,  1989. 

[28]  J.  Sethian,  “Fast  marching  level  set  methods  for  three-dimensional  photolithography  development,”  Proc.  SPIE 
International  Symposium  on  Microlithography,  Santa  Clara,  California,  March,  1996. 

[29]  J.  A.  Sethian,  “A  fast  marching  level-set  method  for  monotonically  advancing  fronts,”  Proc.  Nat.  Acad.  Sci.  93:4, 
pp.  1591-1595, 1996. 

[30]  J.  B.  Tenenbaum,  V.  de  Silva,  and  J.  C.  Langford,  “A  global  geometric  framework  for  nonlinear  dimensionality 
reduction,”  Science,  pp.  2319-2323,  December  2000. 

[31]  R.  Tsai,  L.-T.  Cheng,  S.  Osher,  and  H.-K.  Zhao,  “Fast  sweeping  algorithms  for  a  class  of  Hamilton- Jacobi 
equations,”  UCLA- CAM  Report,  October  2001  (http://www.math.ucla.edu/applied/cam/index.html). 

[32]  J.  N.  Tsitsiklis,  “Efficient  algorithms  for  globally  optimal  trajectories,”  IEEE  Transactions  on  Automatic  Control 
40  pp.  1528-1538, 1995. 

[33]  F.  E.  Wolter,  “Cut  loci  in  bordered  and  unbordered  Riemannian  manifolds,”  Doctoral  Dissertation,  Technische 
Universitat  Berlin,  1985. 

[34]  M.  Zwicker,  M.  Pauly,  O.  Knoll,  and  M.  Gross,  “PointShop  3D:  An  interactive  system  for  point-based  surface 
editing,”  SIGGRAPH,  2002. 


10 


